Zadanie 1¶

Porównać w języku Julia reprezentację bitową liczby 1/3 dla Float16, Float32, Float64 oraz liczby, która jest inicjalizowana jako Float16, a potem rzutowana na Float64.

Rozwiązanie¶

In [17]:
#using Pkg
#Pkg.add("DataFrames")
using DataFrames

fl = [Float16(1/3), Float32(1/3), Float64(1/3), Float64(Float16(1/3))]
fl_bit = map(x -> bitstring(x), fl)

floats = DataFrame()
floats.Float_type=["Float16", "Float32", "Float64", "Float16->Float64"]
floats.Sign_bit = map(x -> x[1], fl_bit)
floats.Exponent_bits = [fl_bit[1][2:6], fl_bit[2][2:9], fl_bit[3][2:12], fl_bit[4][2:12]]
floats.Mantissa_bits = [fl_bit[1][7:16], fl_bit[2][10:32], fl_bit[3][13:64], fl_bit[4][13:64]]

floats
Out[17]:
4×4 DataFrame
RowFloat_typeSign_bitExponent_bitsMantissa_bits
StringCharStringString
1Float160011111111010101
2Float3200111111110101010101010101010101
3Float640011111111010101010101010101010101010101010101010101010101010101
4Float16->Float640011111111010101010101000000000000000000000000000000000000000000

Okazuje się, że po rzutowaniu Float16 na Float64 pozostała część mantysy wypełniana jest zerami. Naturalnie taka liczba mniej dokładnie oddaje wartość 1/3 niż liczba pierwotnie zainicjalizowana jako Float64.

Zadanie 2¶

Zbadać, jak zmienia się odległość między kolejnymi liczbami zminnoprzecinkowymi reprezentowanymi w komputerze za pomocą języka Julia. Narysować wykres używając Plots zależności odległości od wartości liczby dla zakresu od 1.0 do 1000000.0.

Rozwiązanie¶

In [13]:
#using Pkg
#Pkg.add("Plots")
using Plots
r = 1:1000000
epsilons = [eps(Float64(n)) for n in r]
plot(r, epsilons, title="", label="", xlabel="value of the number", ylabel="distance between neighbouring numbers")
Out[13]:

Zadanie 3¶

Jedną z bibliotek numerycznych, jaką dodatkowo będziemy używać na zajęciach jest GSL (język C). Opis jak używać . Korzystając ze wsparcia dla wyświetlania reprezentacji liczb zmiennoprzecinkowych zobaczyć jak zmienia się cecha i mantysa dla coraz mniejszych liczb. Zaobserwować, kiedy matysa przestaje być znormalizowana i dlaczego?

Kod załączyć jako komórka Markdown sformatowana jako C (link). Wynik także jako Markdown (kod albo fragment zrzutu ekranu).

Rozwiązanie¶

Kod w C z użyciem biblioteki gls:

#include <stdio.h>
#include <gsl/gsl_ieee_utils.h>
#include <math.h>

int main()
{
    // floats are 32bit
    float f = 1.12345678 * exp2(-124);

    for (size_t i = 0; i < 27; i++)
    {
    	f/=2;
    	gsl_ieee_printf_float(&f);
    	printf("\n");
    }
    
    return 0;
}

Poniżej wydruk konsoli:

In [14]:
1.00011111100110101101111*2^-125
 1.00011111100110101101111*2^-126
 0.10001111110011010111000*2^-126
 0.01000111111001101011100*2^-126
 0.00100011111100110101110*2^-126
 0.00010001111110011010111*2^-126
 0.00001000111111001101100*2^-126
 0.00000100011111100110110*2^-126
 0.00000010001111110011011*2^-126
 0.00000001000111111001110*2^-126
 0.00000000100011111100111*2^-126
 0.00000000010001111110100*2^-126
 0.00000000001000111111010*2^-126
 0.00000000000100011111101*2^-126
 0.00000000000010001111110*2^-126
 0.00000000000001000111111*2^-126
 0.00000000000000100100000*2^-126
 0.00000000000000010010000*2^-126
 0.00000000000000001001000*2^-126
 0.00000000000000000100100*2^-126
 0.00000000000000000010010*2^-126
 0.00000000000000000001001*2^-126
 0.00000000000000000000100*2^-126
 0.00000000000000000000010*2^-126
 0.00000000000000000000001*2^-126
 0
 0
Out[14]:
0

Mantysa przestaje być znormalizowana, gdy liczba, którą system musi przechować jest mniejsza niż $2^{-126} \cdot 1.0$. Wynika to z tego, że minimalna wartość mantysy to $1$ (dzięki domyślnej jedynce), a najmniejsza wartość cechy bez wszystkich zerowych bitów (aby była znormalizowana) to 1 (00000001), a po uwzględnieniu $biasu=-127$ ostatecznie mamy $2^{1-127}=2^{-126}$. Poniżej tej wartości granicznej liczby muszą być reprezentowane w postaci zdenormalizowanej (cecha ma postać (00000000)) a mantysa traci domyślną jedynkę, przez co należy teraz do przedziału (0, 1). Najmniejszą możliwą do przechowania liczbą w takim układzie jest $(0,00000000, 00000000000000000000001)$

Zadanie 4¶

Na przykładzie wybranego algorytmu niestabilnego numerycznie:

  1. Pokazać, że działa źle.
  2. Pokazać które konkretnie działania powodują zwiększenie błędu (np. dzielenie przez małą liczbę, cancellation).
  3. Zademonstować wersję stabilną.

Wszystkie punkty przedstawić w postaci notatnika Julii.

Rozwiązanie¶

Jako nasze zadanie przyjmijmy obliczenie wartości całek $$y_n = \int_{0}^{1} \frac{x^n}{x+5} dx$$ dla n=0, 1, 2, ... 15.

Zauważmy, że $$ y_n + 5y_{n-1} = \int_0^1 \frac{x^n + 5x^{n-1}}{x+5} dx = \int_0^1 \frac{x^{n-1}(x+5)}{x+5} dx = \left [ \frac{x^n}{n} \right ] _0 ^1 = \frac{1}{n} $$ co daje nam wzór rekurencyjny $$y_n + 5y_{n-1} = \frac{1}{n}$$ Po przekształceniu wzoru otrzymujemy zależność $$y_n = \frac{1}{n} - 5y_{n-1}$$ Potrzebujemy jeszcze wartości pierwszego wyrazu $$y_0 = \int_0^1 \frac{1}{x+5} = \left [ \ln (x+5) \right ] _0^1 = \ln 6 - \ln 5 \approx 0.18232156$$

Przechodzimy do obliczania kolejnych $y_n$

1)¶

In [15]:
numOfSteps = 15
y = Vector{Float32}()
append!(y, 0.18232156)

for n in 1:numOfSteps
    append!(y, 1/n - 5*last(y))
end
results = DataFrame(y_index = [i for i in 0:numOfSteps], value = y)
Out[15]:
16×2 DataFrame
Rowy_indexvalue
Int64Float32
100.182322
210.0883922
320.058039
430.0431383
540.0343086
650.0284571
760.0243812
870.020951
980.020245
1090.00988603
11100.0505698
1211-0.16194
13120.893034
1413-4.38825
151422.0127
1615-109.997

Dla $y_{11}$ otrzymujemy oczywistą sprzeczność, gdyż całka oznaczona w granicach od 0 do 1 z funkcji nieujemnej $\frac{x^n}{x+5}$ nie może mieć wartości mniejszej od zera

2)¶

Powyższy algorytm z każdą iteracją zwiększa błąd, wynikający z niedosknałości reprezentacji liczb zmiennoprzecinkowych w komputerze. Winę za to ponosi wykorzystana przez nas zależność rekurencyjna $$y_n = \frac{1}{n} - 5y_{n-1}$$ Wynika z niej, że z każdą iteracją mnożymy już posiadany błąd razy 5, więc, w szczególności, pierwszy popełniony błąd $\epsilon_0$ (podczas zaokrąglenia $y_0$) będzie już znacząco wpływać na wartość np. wyrażenia $y_8$, bo będzie to $5^8 \cdot \epsilon_0$!

3)¶

Aby algorytm był jednak stabilny, zastosujemy inne podejście. Przekształcamy znaną już zależność rekurencyjną do postaci $$y_{n-1} = \frac{1}{5n} - \frac{1}{5} y_n$$

Teraz błąd będzie z każdym krokiem dzielony przez -5. Ponieważ $y_n$ maleje gdy $n$ rośnie, możemy przypuszczać, że dla dużych $n$, $y_n$ maleje wolno, więc przyjmijmy, że $y_{16} \approx y_{17}$ i korzystając z wzoru $$y_{16} = \frac{1}{5 \cdot 17} - \frac{1}{5}y_{17}$$ otrzymujemy $$y_{16} \approx \frac{1}{5 \cdot 17} - \frac{1}{5} y_{16} \Rightarrow y_{16} \approx \frac{1}{5 \cdot 17} \cdot \frac{5}{6} \approx 0.009803921$$

Teraz stosujemy analogiczne podejście, jak w poprzednim algorytmie, ale zaczynając oczywiście od $y_{16}$ a na $y_0$ kończąc

In [16]:
numOfSteps = 16
y = Vector{Float32}()
append!(y, 0.009803921)

for n in numOfSteps:-1:1
    append!(y, 1/(5*n) - last(y)/5)
end
results = DataFrame(y_index = [i for i in numOfSteps:-1:0], value = y)
Out[16]:
17×2 DataFrame
Rowy_indexvalue
Int64Float32
1160.00980392
2150.0105392
3140.0112255
4130.0120406
5120.0129765
6110.0140714
7100.0153675
890.0169265
980.0188369
1070.0212326
1160.0243249
1250.0284684
1340.0343063
1430.0431387
1520.0580389
1610.0883922
1700.182322

Tym razem wartość $y_0 \approx 0.182322$, co jest poprawnym wynikiem

In [ ]: